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Abstract 


A three-dimensional multiblock Navier-Stokes code, PAB3D, which was 
developed for propulsion integration and general aerodynamic analysis, has 
been used extensively by NASA Langley and other organizations to perform both 
internal (exhaust) and external flow analysis of complex aircraft configurations. 
This code was designed to solve the simplified Reynolds Averaged Navier- 
Stokes equations. A two-equation k-e turbulence model has been used with 
considerable success, especially for attached flows. Accurate predicting of 
transonic shock wave location and pressure recovery in separated flow regions 
ha s been more difficult. Two algebraic Reynolds stress models (ASM) have 
been recently implemented in the code that greatly improved the code's ability 
to predict these difficult flow conditions. Good agreement with Direct 
Numerical Simulation (DNS) for a subsonic flat plate was achieved with ASMs 
developed by Shih, Zhu, and Lumley and Gatski and Speziale. Good predictions 
were also achieved at subsonic and transonic Mach numbers for shock location 
and tr ailin g edge boattail pressure recovery on a single-engine afterbody /nozzle 
model. 



Introduction 


CFD methods along with accurate turbulence models are required to predict 
aerodynamic effects at transonic conditions. Accurate prediction of pressure 
distribution and skin friction coefficient is of paramount importance to the design 
of aerodynamic configurations. Accurate predictions of boundary layer structure 
and flow separation by CFD methods are also very critical. 

It is widely accepted that the computational economy of two-equation 
turbulence models (ref. 1) offers a reasonable compromise for computing practical 
flow problems. The theoretical advantages of these models over algebraic models 
are the incorporation of turbulence history-dependent non-local effect (through 
the convection and viscous diffusion of the Reynolds stress) which are known to 
play an important role in determining the turbulence structure in complex flows. 
In the standard two-equation k-e turbulence model, transport equations are 
carried for kinetic energy and dissipation (ref. 1). The k-e equations can be 
applied to the near wall region as well as far away from wall boundaries. For 
regions of the flow far away from solid boundaries, the high Reynolds number form 
of the model can be used while wall damping functions must be used when 
applying the model near wall boundaries. 

Reynolds Stress models (refs. 2 & 3) have the potential of producing more 
accurate turbulent simulation as compared with two-equation turbulence models. 
However, numerical calculations using the more advanced Reynolds Stress models 
require the solution of transport equations for each individual component of the 
Reynolds stress tensor (five more equations) besides solving the Navier-Stokes 
equations. This approach requires tremendous computational time when solving 
three-dimensional (3-D) flow problems. Moreover, the transport equations for the 
second-order models require closure approximation for higher-order turbulence 
correlation that has uncertain physical foundations. 

As an alternative to Reynolds Stress turbulence models, Shih et al. 4 and 
Gatski and Speziale 5 developed algebraic forms of the Reynolds stress model, 
which broadened the range of applicability of the existing linear models (k-e for 
example) while maintaining most of their popular features (such as reduction to 
mixing layer theories for thin shear flows and the possibility of implementation in 
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existing Navier-Stokes codes without substantially increasing the computational 
time). These algebraic stress models (ASM) were developed by making an 
asymptotic expansion subject to constraints of dimensional and tensorial 
invariance, realizability, and material frame indifference. The resulting ASM 
models showed substantially improved predictions in incompressible turbulent 
channel flows and yielded normal Reynolds-stress differences that gave rise to 
secondary flows in complex flows. 

In the present research work, the ASM models developed by Shih, Zhu, and 
Lumley (SZL) and Gatski and Speziale (GS) were implemented and tested in the 
CFD code PAB3D. The GS algebraic stress model uses coefficients derived from 
the Sarker, Gastki, and Speziale^ (SSG) Reynolds Stress model. PAB3D is a 
general purpose, 3-D, multiblock Navier-Stokes code described and applied in 
ref. 1. The flow solver contains the Baldwin-Lomax 6 turbulence model and a two- 
equation k-e turbulence model with various near wall damping functions. The two 
ASM models identified above were implemented for the present study to resolve 
flow field anisotropy. The code has a built in performance module to compute 
various quantities such as lift, drag, thrust, and discharge coefficient. During a 
typical numerical simulation, these quantities are constantly monitored to assess 
the performance of the propulsion system under consideration. 

Computed results using ASM are compared with published experimental 
data at several levels of increasing flow complexity. The geometries considered in 
the present study are a subsonic flat plate 7 , a transonic axisymmetric afterbody 8 , 
and a supersonic square duct 9 . The flat plate was selected because it is the 
simplest of all the geometries for which experimental and direct numerical 
simulation 10 (DNS) data are available. The axisymmetric afterbody represents the 
next level of flow complexity because this case contains a separated flow region. 
The separation is generated by either an adverse pressure gradient or a shock 
boundary layer interaction. The square duct case provides a simple, full 3-D test 
case with complex flow structure. In this case, secondary flow is generated by the 
anisotropy of the normal stress components. 
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Sy mb ols and Abbreviations 


C r C 2 , Cg.and C 4 
C el’ and C e2 

Cn 

Cp 

D 

k 

1 

n 

Sy 

Uy 


Rex 


u.v.w 

u+ 

x,y,z 

8 

K 

e 

Y 

P 

v 

P 

Xw 

Q 

(t ii. 0 


model constants for the GS Algebraic Stress model 
model constants for the k-e turbulence model 

turbulence viscosity coefficient for the k-e model . .09 

pressure coefficient 

duct width, cm. 

turbulent kinetic energy. Pa. 

length of the axisymmetric model 

normal distance from the wall, cm. 

Strain component, 1/sec 

— i-, velocity gradient, 1/sec 
dxj 

law-of-the-wall coordinate , n A/Pw T w _ n Pw u _x 

M-w Pw 

Reynolds number based on streamwise distance from 

the plate leading edge, u^p^/Poo- 

velocity components in the (x t y,z) directions, m/sec. 

u _ Re c 

law-of-the wall coordinate, — — r - 

u T n 

friction velocity, ^t w /p 
spatial coordinates, cm. 

Kronecker delta and boundary layer thickness 

von Karman constant. 

turbulent energy dissipation 

ratio of specific heat, 1.4. 

dynamic viscosity coefficient, m^/sec. 

kinematic viscosity coefficient, m 2 /sec. 

density, kg/m 3 

Reynolds stress components. 

9u 

wall shear stress, |i— l w 

dn 

vorticity 

generalized coordinates as functions of (x,y,z) 


Superscripts 

L laminar 

X turbulent 
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Abbreviations 


ASM algebraic stress model 

DNS direct numerical simulation 

GS Gatski and Speziale algebraic stress model 

SSG Sarker, Gatski, and Speziale Reynolds stress model 

SZL Shih, Zhu and Limley algebraic stress model 

k-e k-epsilon model - Jones-Launder wall damping 

V13 Version 13 of PAB3D code 


Computational Procedure 


Code Structure 


Figure 1 presents a sketch of the logical structure of the PAB3D code. This 
code solves the Reynolds Averaged Navier-Stokes (RANS) equations, along with 
several different forms of a two-equation turbulence model, written in a general 
3-D form. It requires 15 words per grid point to compute laminar or turbulent 
flows using algebraic turbulence models and 20 words per grid point to compute 
turbulent flows using a two-equation turbulence model. The code speed is 33jis 
and 43ps per grid point (Cray 2 time) for solving laminar and turbulent flows, 
respectively. The convective (inviscid) terms of the transport equations (RANS 
and turbulence equations) are simulated using state-of-the-art upwind schemes 1. 
The viscous and diffusion terms in these equations are evaluated using a central 
difference scheme. The two-equation turbulence model equations are uncoupled 
from the RANS equations and solved with a different time step. Different near- 
wall, compressibility and damping coefficients can be used with any of the two- 
equation turbulence model options available in the code. Recently, several 
different forms of Algebraic Reynolds Stress turbulence models (ASM) 4 ’ 5 have 
been implemented in the PAB3D code (VI 3). The code utilizes post- and pre- 
processing software packages for manipulation of the computational grids, 
solutions, and the creation of the patched block interface data base. A 3-D 
performance package is utilized to calculate aerodynamic quantities such as nozzle 
discharge coefficient and thrust ratio for exhaust systems, and lift and drag 
coefficient for aerodynamic bodies. 
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Eddv Viscosity Turbulence Models 

In most 3-D Navier-Stokes codes, the stress components are modeled as: 
T = X- -I- 1- 

W IJ fc lJ ' W 1J 

Where, T is calculated with the eddy viscosity formulation through the 
Boussinesq's hypothesis and following Kolmogorov: 
xJ = A T 8 y -2mSy 


A T =-(pk + p x S kk ) 

A L =fn L s tt 

k 2 

Mr - — 

£ 

and, fjj. is an empirical function and Cjj. is set to a value of .09 for the standard k-e 
model (linear model). K is the turbulent kinetic energy and e is the rate of 
dissipation. Despite the current popularity of the k-e turbulence model, it still 
shares a common deficiency with the mixing length algebraic turbulence models, 
i.e. the Boussinesq eddy viscosity hypothesis. It is well known that these models 
give inaccurate prediction of normal Reynolds stress differences, which makes the 
eddy viscosity models incapable of simulating complex flows. In the above 
formulation, stress components are directly related to the rate of strain S after 
dropping the k term. This simplifies the implementation of the two equation 
turbulence model in most 3-D Navier-Stokes codes due to the eddy viscosity form 
as described in Ref. 1 . 


gebraic Stress Models 


The use of the algebraic stress models (ASMs) is not as straight forward as 
the implementation of lower level turbulence models (including the two-equation 
linear turbulence model). Because of the nonlinear relation between stress and 
the velocity gradients, turbulence stresses must be evaluated at the cell faces 
instead of the cell centers. The six stress components (laminar and turbulent) are 
computed for each cell. The full 3-D codes require the calculation of each stress 
component at the six faces of each cell for a total of 36 values. Fortunately, in 
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most cases, the viscous diffusion terms can only be evaluated in one or two 
directions, mostly due to grid restrictions. This reduces the number of variables 
to 12 or 24. 

In the present study, two ASM models (SZL and GS) were implemented in 
the PAB3D code. SZL 4 developed a new Reynolds stress algebraic equation model 
using a general turbulent constitutive relation (Shih and Lumley 11 ). In this 
development, the constraints based on rapid distortion theory and realizability are 
imposed. In this model, turbulent stress components are composed of linear and 
nonlinear parts as follows: 


t T _ _T(1) , T(n) 
l ij L ij ' L ij 

xf ) = A5 ij -2n,S ij 
xX ( ”»=2|}^(Q* a[ s y-S; k Q; j ) 


where, 

P J 


'l-9cj(— ) 2 

H e 


1+6 


S*k£2*k 


Cji = 


e e 
1 


6.5+A 


*U*k 


A* = V6 cos(<))) 

<|> =icos _1 (V6W*) 

^ij — ^ij — 2 ^ U i-j ~ U j'> ^ = _ ^ji 


sXs* 


w ‘ = and U* = ^s; + Q*yQ* y 


GS 5 describes the relationship between the new explicit algebraic stress 
models and anisotropic eddy viscosity models in reference 5. The GS form of the 
explicit algebraic stress model was implemented in PAB3D as follows: 
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+ p2(S ik S tj -iS mI ,S mn 8 ij )] 

p _ 3(l+ti 2 )a 

^ 3+ti 2 +6ti' z ; 2 +6; 2 

p,=g(2-c 4 ) and p 2 = 2g(2-c 3 ) 

«-g(j-c 2 ) 

n 2 =s; j s; j ,e = n: j n- 

S ‘ J =^g( 2 - C 3)S U 

n: j = - g ( 2 -c 4 )n ij 

J 2£ J 

1 

Cl , ^2-1 
2 C el -1 

C, = 3.4 + 5^1, C 2 = 0.368, C 3 = 1.25,C 4 =0.4 

C El -1 

Results and Discussion 

Because of the complex nature of the flow discussed in this paper, it was 
essential to maintain appropriate grid spacing near solid boundaries to ensure 
appropriate near wall behavior of k and e. For the current study, y + of the first 
point located off the wall was less than one and the grid was stretched in the 
normal direction using an exponential grid stretching formula. Approximately 16 
points were placed normal to the walls to resolve the boundary layer. The rest of 
the points in the normal direction were distributed uniformly between the edge of 
the boundary layer and the farfield or the symmetry boundary. 

Subsonic Flow Over a Flat Plate 

The first test case considered in the present study was subsonic flow over a 
flat plate (ref. 7). Four physical properties of flow over a flat plate were chosen to 
evaluate the accuracy of the different forms of the two-equation and algebraic 
stress models. These parameters were evaluated at Re of 1420. The four 

properties are: 

1) u + (normalized velocity) 
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+ * / 

U = U / U T 

Uoo 

U T = J Vp / Pwdu 
0 

This parameter is a standard variable in evaluating turbulence model accuracy. 
Laminar sublayer, buffer, and turbulent regions should be clearly defined during 
this evaluation. 

2) k + (normalized turbulence kinetic energy) 

k + = k/u? 

3) uV’ + (normalized shear stress) 
u' v ,+ = u' V /u 2 = X 12 t u 2 

4) u ,+2 , v’ +2 , w' +2 (normalized normal stress components), for example 

u ,+2 =u ,2 /u 2 = T n /u 2 

These parameters were selected to show three-dimensionality of normal stress 
components for a two-dimensional (2-D) flow case. 

Figure 2 shows the comparison of the calculated boundary layer parameters 
with DNS 10 and experimental data. In figure 2a, the standard Jones and Launder 
two-equation turbulence model shows substantial deficiencies in predicting any of 
the four parameters. In predicting the normalized velocity, u + , the flow 
transitioned from a laminar sublayer directly to a fully turbulent region without 
going through a buffer region. The model under predicted the k+ peak by about 
20 percent. Also, the peak shear stress was under predicted by at least 10 
percent. As expected, the turbulence model did not predict the three- 
dimensionality of the normal stress components due to the two-equation k-e 
turbulence model assumptions of a linear relation between stress and strain. 

The second turbulence model evaluated is the SZL algebraic stress model. 
Figure 2b shows the evaluation of this turbulence model in prediction of the four 
properties. The u+ is in excellent agreement with experimental as well as the 
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theoretical data. All of the flow regions are well defined in this prediction. The 
model accurately predicted the k+ peak as compared with the DNS data. 

Similarly, the shear stress peak is more accurately predicted although somewhat 
low. The model also shows the anisotropy of the three normal stress components. 

The algebraic stress model based on the SSG coefficients are compared with 
the data in figure 2c. This model provides the best agreement with DNS and 
experimental data when compared with the other two turbulence models. Correct 
kinetic energy and shear stress peaks and the general features of all four 
parameters are captured well. However, this turbulence model required a smaller 
time step to avoid divergence of the solution and it needed 20 percent more 
iterations to achieve similar solution accuracy as the SZL model. 

Transonic Axisymmetric Afterbody /Nozzle 

The second test case considered was an axisymmetric afterbody/ nozzle with 
a constant diameter cylinder representing a solid jet exhaust plume simulation 
(see ref. 8). This configuration, although axisymmetric, can involve solutions 
where all three dimension components of Reynolds stress are non-zero and 
unequal. Mach numbers of 0.6, 0.85 and 0.9 were the selected freestream flow 
conditions for this test case. For this configuration, a separated flow region occurs 
on the nozzle boattail through different flow mechanisms depending on 
freestream Mach number. At M = 0.60, separation is caused by an adverse 
pressure gradient. The separated flow regions at Mach numbers of 0.85 and 0.90 
occur due to a shock/boundary-layer interaction. 

Figure 3 shows the comparison between the experimental data (ref. 8) for 
this test case and two different turbulence model predictions. The standard k-e 
model and the GS form of the algebraic Reynolds stress model are shown for this 
prediction. There was little difference between results using the SZL (not shown) 
or GS models except for convergence history. All the calculations agree well with 
the data up to the separated flow region or the location of the shock. In general, 
the standard (linear) two-equation turbulence model overpredicted the pressure 
recovery on the afterbody by at least 40 percent as shown in figure 3a. For the 
Mach = 0.90 case, the linear turbulence model predicted the shock location 
downstream of the data. The GS form of the algebraic Reynolds stress (nonlinear) 
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model accurately predicted the pressure recovery in the separated region as well 
as the calculated shock location (see figure 3b). 

Supersonic Square Duct 

The final test case geometry was a square duct (ref. 9). Numerical 
calculations were carried out for supersonic flow through square duct using linear 
and SZL k-e turbulence models. The duct Mach number was 3.9 and the unit Re 
was 1.2 million/m. Figure 4 shows the cross flow velocity patterns computed 
using the two different turbulence models at X = 20, 30, 40 and 50. Dramatically 
improved results are obtained using the SZL model as shown in figure 4b. Similar 
to the experimental results reported in reference 9, the SZL computations clearly 
show that the secondary flows (vortices) are symmetrical about the diagonal and 
rotate in opposite directions. These vortices are essentially driven by the 
gradients of the Reynolds stresses which transport net momentum towards the 
comer of the duct and cannot be simulated by the linear models. The computed 
cross flow velocity vectors using the SZL turbulence model are in good agreement 
with the experimentally observed patterns. In contrast, the linear model predicts 
a unidirectional flow due to the inability of the turbulence model to represent the 
flow physics. 

Concluding Remarks 

An investigation of the effect of using standard k-e and algebraic Reynolds 
stress turbulence models on the computed flow field of several different 
aerodynamic configurations was conducted. The computed results were 
compared with the available computations, experimental, and theoretical data 
bases. The geometries considered in the present study were a flat plate, a 
transonic axisymmetric afterbody and supersonic flow through a square duct. 

The flat plate test case was selected because it is the simplest of all the 2-D 
geometries which show the three-dimensionality of the normal turbulent stress 
components. The results obtained using either the SZL or the GS form of the 
explicit algebraic Reynolds stress (ASM) turbulence model were more accurate 
and consistent than the standard (linear) k-e turbulence model. The ASM 
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turbulence model accurately predicted the turbulent kinetic energy and shear 
stress peaks as well as the anisotropy of the normal stress components. 

The axisymmetric afterbody geometry represented the next level flow 
complexity because this case contains separated flow regions driven by either an 
adverse pressure gradient or a shock/boundaiy-layer interaction. The algebraic 
Reynolds stress turbulence model generally provided accurate predictions of the 
pressure recovery in separated flow regions whereas the standard k-e turbulence 
model overpredicted the pressure recovery in these regions. In addition, the 
ASM turbulence model generally provided better predictions of shock location. 

For the test case of square duct geometry, a secondary flow structure 
developed in directions perpendicular to the main flow. The ASM results clearly 
showed that the secondary flows (vortices) were symmetrical about the duct 
diagonal and rotate in opposite directions. These vortices are essentially driven by 
the gradients of the Reynolds stresses, which could not be simulated by the linear 
model. 

This investigation provided significant insight into the applications of 
explicit algebraic Reynolds stress turbulence models in predicting attached and 
separated flows. The results of this investigation support a high level of 
confidence in using these advanced turbulence models for more complex 
aerodynamic configurations. 
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Figure 1. Sketch of Logical Structure of PAB3D code. 








(a) Standard Jones & Launder 2-eqn. k-e 

Figure 2.- Comparison of calculated boundary layer parameters with DNS and data. 
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(b) Shih, Zhu & Lumley algebraic stress model 
Figure 2.- Continued. 
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(a) standard Jones & Launder 2-eqn. k-e 

Figure 3. Experimental and predicted pressure distribution on axisymmetric 
afterbody with solid plume. 
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(b) Gatski & Speziale algebraic stress model. 
Figure 3. Concluded 


o 

CO 

II 

X 


IT) 



L-L “ ‘ 11 ■ iLAliJUJ-UJ iimuHim 

m ^ q oi r- cp 

d o d o o o 


Q 

7a 


o 

in 

ll 

X 


in 
d 

"■ >» ^ ^ ^ 

>* N N ^ X \ 

* ' N N S. \\\XV^^^^^^^^^^^^i £T) 

- ' ^ ^ N * 

\ \ N \ \ N N\NNX\N^^^\\\\a o 

n \ \ \ \. 

\ \ \ \ \ \ 

\ ^ \ \\\\\n 

\ \ \ \\\\\> 

° 

ii mu ° 




0 

> 


d 

c\i 

il 

X 


in 



d> d> ci d> cd d> 


Q 

7a 


o 

^r 

n 

X 



.^^^^^^\\\\\\\\\^* 
v\\\\\\\\\\\\\W 

X\\\\\\\\\\\w 
xwswwww 
,\\\\\\\\u« 

WWNWWW 
WWWWWW 
. AWWWWW 
. AWWWUWWI 
AWWWVWWW 
, , AWNNSNWVttrt 
J\S\N\\\S\\U\\W 
\\W\\\\\WN\\«r 


m 

d 

d 

co 

d 

CM 

d 


o 

o 


Q 

> 


21 


(a) Standard Jones & Launder 2-eqn. k-e. 

Figure 4. Cross-flow velocity pattern for supersonic square duct. 



CO 


(b) Shih, Zhu & Lumley algebraic stress model 
Figure 4. Concluded 
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